Sesión 04

Diseño Cuadrado Latino


Christian Vásquez-Velasco, Bach., M.Sc.(c)

InkaStats Academy

2023

Planeamiento


Instalar paquetes necesarios


# install.packages("devtools")
# devtools::install_github("emitanaka/edibble", force = T)
# devtools::install_github("emitanaka/deggust", force = T)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(readxl, agricolae, agricolaeplotr, car, tidyverse, PMCMRplus, outliers, nortest, mvtnorm, lmtest, ExpDes, edibble, gt,
               gtsummary, devtools, deggust, xlsx, desplot,
               ggResidpanel, fastGraph, gvlma,
               ScottKnott, kSamples, psych)

Crear un libro con el paquete agricolae


variedades <- c("Amazon", "Coari x Lame", "Coari x Yangambí",
                "Unipalma")

salida <- agricolae::design.lsd(trt = variedades,
                                serie = 2,
                                seed = 123,
                                kinds = "Mersenne-Twister",
                                first = TRUE,
                                randomization = TRUE)

Sketch del diseño

salida$sketch
     [,1]               [,2]               [,3]              
[1,] "Coari x Lame"     "Coari x Yangambí" "Unipalma"        
[2,] "Amazon"           "Coari x Lame"     "Coari x Yangambí"
[3,] "Coari x Yangambí" "Unipalma"         "Amazon"          
[4,] "Unipalma"         "Amazon"           "Coari x Lame"    
     [,4]              
[1,] "Amazon"          
[2,] "Unipalma"        
[3,] "Coari x Lame"    
[4,] "Coari x Yangambí"
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Sketch codificado del diseño

print(matrix(salida$book[,1], byrow=TRUE, ncol=4))
     [,1] [,2] [,3] [,4]
[1,]  101  102  103  104
[2,]  201  202  203  204
[3,]  301  302  303  304
[4,]  401  402  403  404

Guardar el libro generado

write.table(salida$book, "books/lsd.txt",
            row.names = FALSE, sep = "\t")

write.xlsx(salida$book, "books/lsd.xlsx", sheetName = "book", append = FALSE, row.names = FALSE)

Libro de campo

fieldbook <- salida %>% 
  zigzag()

fieldbook %>%
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro de campo generado

write.table(fieldbook,
  "books/lsd.txt",
  row.names = FALSE, sep = "\t")

write.xlsx(fieldbook,
  "books/lsd.xlsx",
  sheetName = "book",
  append = FALSE,
  row.names = FALSE)
agricolaeplotr::plot_latin_square(design = salida,
                                  factor_name = "variedades",              
                                  reverse_y = TRUE,
                                  reverse_x = FALSE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas",
       title = "Diseño Cuadrado Latino",
       subtitle = "Rendimiento de cultivos",
       caption = "Curso: Métodos Estadísticos para Experimentos Agrícolas en R - 2023") +
  theme_classic()

Crear un libro de campo con el paquete edibble


menu_lsd()
design("Latin Square Design") %>%
  set_units(row = 5,
            col = 5,
            unit = crossed_by(row, col)) %>%
  set_trts(trt = 5) %>%
  allot_trts(trt ~ unit) %>%
  assign_trts("random", seed = 590) %>%
  serve_table()
lsd <- takeout(menu_lsd(t = 4,
                        seed = 123))
lsd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
lsd2 <- design("Latin Square Design") %>%
  set_units(row = 5,
            col = 5,
            unit = crossed_by(row, col)) %>%
  set_trts(variedades = c("Amazon", "Coari x Lame", "Coari x Yangambí", "Unipalma")) %>%
  allot_trts(variedades ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
lsd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
deggust::autoplot(lsd)

Análisis del diseño


Importación de datos


archivos <- list.files(pattern = "ejercicio 1.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor)

Creación del modelo lineal


modelo.dcl <- lm(rdt ~ trt + `F` + `C`, data = data)

Definición del modelo


\[Y = \beta_0 + \beta_1*T_B + \beta_2*T_C + \beta_3*T_D + \beta_4*T_E + \beta_5*F{2} + \beta_6*F{3} + \beta_7*F{4} + \beta_8*F{5} + \beta_9*C{2} + \beta_{10}*C{3} + \beta_{11}*C{4} + \beta_{12}*C{5} + \epsilon\]

\[\hat{Y} = \beta_0 + \beta_1*T_B + \beta_2*T_C + \beta_3*T_D + \beta_4*T_E + \beta_5*F{2} + \beta_6*F{3} + \beta_7*F{4} + \beta_8*F{5} + \beta_9*C{2} + \beta_{10}*C{3} + \beta_{11}*C{4} + \beta_{12}*C{5}\]

summary(modelo.dcl)

Call:
lm(formula = rdt ~ trt + F + C, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9796 -0.9756  0.1684  0.7484  3.1784 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   44.670      1.567  28.515 2.15e-12 ***
trtB           4.584      1.374   3.336  0.00593 ** 
trtC           8.594      1.374   6.255 4.22e-05 ***
trtD          12.034      1.374   8.759 1.47e-06 ***
trtE          11.744      1.374   8.548 1.90e-06 ***
FF2            2.152      1.374   1.566  0.14326    
FF3           -1.864      1.374  -1.357  0.19985    
FF4           -1.576      1.374  -1.147  0.27371    
FF5            3.154      1.374   2.296  0.04052 *  
CC2            1.862      1.374   1.355  0.20031    
CC3            2.006      1.374   1.460  0.16996    
CC4            2.864      1.374   2.085  0.05915 .  
CC5            3.724      1.374   2.710  0.01894 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.172 on 12 degrees of freedom
Multiple R-squared:  0.921, Adjusted R-squared:  0.8419 
F-statistic: 11.65 on 12 and 12 DF,  p-value: 7.933e-05

Conclusión: Se concluye que los coeficientes \(\beta_0\), \(\beta_1\), \(\beta_2\), \(\beta_3\), \(\beta_4\), \(\beta_8\), \(\beta_12\) correspondientes a los tratamientos B, C, D y E, a la fila 5 y columna 5 respectivamente, fueron estadísticamente distintos de 0. El ajuste del modelo descrito en el \(R^2\) fue de 0.921, lo que quiere decir que el 92.1 % de la varianza del rendimiento está explicado por las variables Tratamiento, Fila y Columna.

Verificación visual de los supuestos del modelo


performance::check_model(modelo.dcl)
ggResidpanel::resid_panel(modelo.dcl)
influence.measures(modelo.dcl)
Influence measures of
     lm(formula = rdt ~ trt + F + C, data = data) :

    dfb.1_  dfb.trtB  dfb.trtC  dfb.trtD  dfb.trtE   dfb.FF2   dfb.FF3
1   1.5219 -1.07e-16  3.48e-16  5.77e-16  1.08e+00 -1.08e+00 -1.08e+00
2  -0.0391 -7.48e-17 -7.43e-02 -4.18e-18 -2.87e-17 -7.43e-02 -7.26e-18
3  -0.5548 -1.05e+00  7.63e-16  2.26e-16 -1.57e-16  6.22e-16 -1.05e+00
4  -0.3441  2.45e-01  2.45e-01  2.45e-01  2.45e-01  7.77e-17  1.41e-16
5   0.1609  1.61e-17  2.19e-17  3.06e-01 -3.62e-17  3.15e-17 -4.38e-17
6   0.3081 -2.20e-01 -2.20e-01 -2.20e-01 -2.20e-01 -2.20e-01 -2.20e-01
7  -0.0176  3.06e-17  1.79e-17  5.01e-02  5.49e-18  5.01e-02  2.82e-17
8  -0.0469  4.17e-17  1.34e-01  2.45e-17  9.38e-18 -2.32e-17  1.34e-01
9  -0.0616  1.75e-01 -5.68e-17 -8.73e-18 -1.08e-17 -1.20e-17  2.52e-17
10  0.2180  3.27e-17  4.45e-17  0.00e+00 -6.21e-01  6.23e-17 -2.61e-17
11 -0.1284  8.46e-18 -2.44e-01 -5.90e-17 -1.05e-16  2.44e-01  2.44e-01
12 -0.4084  1.16e+00 -4.90e-16  3.30e-17 -1.02e-16  1.16e+00  1.52e-18
13  0.1012 -1.26e-16 -1.45e-16 -2.88e-01 -7.72e-17 -8.00e-17 -2.88e-01
14  0.1009  2.46e-17  1.46e-17  4.37e-17 -2.88e-01 -1.51e-17 -2.06e-17
15 -0.0612  1.16e-01  1.16e-01  1.16e-01  1.16e-01 -6.53e-17 -5.40e-17
16 -0.3266 -6.21e-01 -1.85e-16 -2.07e-16 -2.67e-16  6.21e-01  6.21e-01
17  0.1322 -9.79e-17  2.70e-17  0.00e+00 -3.77e-01 -3.77e-01  3.76e-17
18  0.3959 -7.52e-01 -7.52e-01 -7.52e-01 -7.52e-01  1.57e-16  7.52e-01
19 -0.0770  4.65e-17  9.44e-17  2.20e-01  2.99e-17  4.23e-17  4.72e-17
20 -0.0172  6.79e-18  4.89e-02  1.39e-18 -4.29e-18 -6.88e-18 -1.03e-18
21 -0.1521 -3.76e-17 -3.11e-16 -2.89e-01 -1.46e-16  2.89e-01  2.89e-01
22 -0.2997  5.69e-01  5.69e-01  5.69e-01  5.69e-01 -5.69e-01 -5.41e-17
23 -0.1230  6.35e-17  7.69e-17  8.00e-17  3.51e-01 -6.08e-17  3.51e-01
24 -0.0469  1.85e-17  1.34e-01 -2.70e-17 -2.46e-17 -2.80e-17 -1.92e-17
25 -0.1232  3.51e-01 -7.86e-17  4.39e-18 -7.40e-17  1.46e-16  2.22e-16
     dfb.FF4   dfb.FF5   dfb.CC2   dfb.CC3   dfb.CC4   dfb.CC5  dffit  cov.r
1  -1.08e+00 -1.08e+00 -1.08e+00 -1.08e+00 -1.08e+00 -1.08e+00  2.473 0.0296
2   2.52e-18 -3.26e-19  7.43e-02  7.43e-02  7.43e-02  7.43e-02 -0.169 6.2578
3   3.18e-16  1.85e-17  1.05e+00  1.05e+00  1.05e+00  1.05e+00 -2.404 0.0378
4  -2.45e-01  7.96e-17  2.45e-01  2.45e-01  2.45e-01  2.45e-01 -0.559 4.6111
5   0.00e+00  3.06e-01 -3.06e-01 -3.06e-01 -3.06e-01 -3.06e-01  0.697 3.8386
6  -2.20e-01 -2.20e-01  2.20e-01 -1.76e-16 -1.41e-16 -1.40e-16  0.501 4.9256
7   2.21e-17  2.37e-17  5.01e-02  1.23e-17  9.26e-18  3.30e-18  0.114 6.3655
8   9.85e-18  1.06e-17  1.34e-01  5.47e-17  6.20e-17  7.33e-17  0.305 5.8371
9   1.75e-01 -3.23e-17  1.75e-01  8.16e-17  9.33e-17  1.23e-16  0.400 5.4280
10  1.89e-16 -6.21e-01 -6.21e-01  2.23e-17 -1.34e-18  6.82e-17 -1.417 0.8533
11  2.44e-01  2.44e-01 -1.27e-16 -2.44e-01 -2.41e-17 -4.82e-17 -0.556 4.6262
12  1.58e-17 -6.13e-17 -1.21e-16  1.16e+00  3.57e-17 -1.28e-16  2.655 0.0154
13 -9.85e-18 -7.34e-17 -1.20e-16 -2.88e-01 -1.91e-16 -1.96e-16 -0.658 4.0624
14 -2.88e-01  2.53e-18  4.99e-17 -2.88e-01 -4.23e-17 -1.89e-17 -0.656 4.0702
15 -6.19e-17 -1.16e-01  1.21e-17 -1.16e-01  1.88e-17  3.32e-17 -0.265 5.9806
16  6.21e-01  6.21e-01 -1.29e-16  0.00e+00 -6.21e-01  0.00e+00 -1.415 0.8570
17  8.83e-17  8.93e-17  6.54e-17  0.00e+00 -3.77e-01  4.14e-17 -0.859 2.9544
18  2.03e-16  1.52e-16 -1.96e-16 -5.39e-17  7.52e-01 -8.25e-17  1.716 0.3662
19  2.20e-01  4.43e-17  8.85e-17  1.10e-16  2.20e-01  9.15e-17  0.501 4.9256
20  7.44e-18  4.89e-02  1.28e-17  1.05e-17  4.89e-02  1.07e-18  0.112 6.3697
21  2.89e-01  2.89e-01 -1.65e-16 -9.14e-17 -8.79e-17 -2.89e-01 -0.659 4.0545
22 -7.22e-17  2.00e-17 -2.09e-17 -1.20e-17  8.66e-17 -5.69e-01 -1.299 1.1557
23 -7.65e-17 -2.77e-17 -1.46e-16  0.00e+00  0.00e+00  3.51e-01  0.800 3.2734
24  1.34e-01  1.29e-17 -1.86e-17  1.63e-17  2.03e-17  1.34e-01  0.305 5.8371
25  2.14e-16  3.51e-01  8.62e-18  5.19e-17  8.01e-17  3.51e-01  0.801 3.2656
    cook.d  hat inf
1  0.33917 0.52   *
2  0.00240 0.52   *
3  0.32660 0.52   *
4  0.02556 0.52   *
5  0.03921 0.52    
6  0.02061 0.52   *
7  0.00109 0.52   *
8  0.00773 0.52   *
9  0.01326 0.52   *
10 0.14417 0.52    
11 0.02532 0.52   *
12 0.37164 0.52   *
13 0.03501 0.52    
14 0.03487 0.52    
15 0.00587 0.52   *
16 0.14387 0.52    
17 0.05837 0.52    
18 0.19808 0.52    
19 0.02061 0.52   *
20 0.00104 0.52   *
21 0.03516 0.52    
22 0.12395 0.52    
23 0.05091 0.52    
24 0.00773 0.52   *
25 0.05108 0.52    
influenceIndexPlot(modelo.dcl)

Cumplimiento de supuestos del modelo lineal general


Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dcl,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1     0.029851793      1.752972  0.2000
   2    -0.456648318      2.721108  0.1656
   3    -0.125539304      1.877683  0.9476
   4    -0.123558763      1.801953  0.7992
   5     0.001742879      1.515653  0.2376
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dcl, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dcl
DW = 1.753, p-value = 0.2041
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dcl))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dcl)
W = 0.96303, p-value = 0.4782
nortest::ad.test(rstudent(modelo.dcl))

    Anderson-Darling normality test

data:  rstudent(modelo.dcl)
A = 0.38392, p-value = 0.3688
lillie.test(rstudent(modelo.dcl))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dcl)
D = 0.13219, p-value = 0.3133
ks.test(rstudent(modelo.dcl), "pnorm",
        alternative = "two.sided")

    Exact one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dcl)
D = 0.10352, p-value = 0.9264
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dcl))

    Cramer-von Mises normality test

data:  rstudent(modelo.dcl)
W = 0.054267, p-value = 0.4357
pearson.test(rstudent(modelo.dcl))

    Pearson chi-square normality test

data:  rstudent(modelo.dcl)
P = 7.96, p-value = 0.1585
sf.test(rstudent(modelo.dcl))

    Shapiro-Francia normality test

data:  rstudent(modelo.dcl)
W = 0.95883, p-value = 0.3324

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento.

ncvTest(modelo.dcl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.915634, Df = 1, p = 0.33862
bptest(modelo.dcl)

    studentized Breusch-Pagan test

data:  modelo.dcl
BP = 15.008, df = 12, p-value = 0.241
bptest(modelo.dcl, studentize = F)

    Breusch-Pagan test

data:  modelo.dcl
BP = 12.906, df = 12, p-value = 0.3759
olsrr::ols_test_breusch_pagan(modelo.dcl)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

             Data               
 -------------------------------
 Response : rdt 
 Variables: fitted values of rdt 

       Test Summary         
 ---------------------------
 DF            =    1 
 Chi2          =    0.915634 
 Prob > Chi2   =    0.338624 

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dcl %>% gvlma()

Call:
lm(formula = rdt ~ trt + F + C, data = data)

Coefficients:
(Intercept)         trtB         trtC         trtD         trtE          FF2  
     44.670        4.584        8.594       12.034       11.744        2.152  
        FF3          FF4          FF5          CC2          CC3          CC4  
     -1.864       -1.576        3.154        1.862        2.006        2.864  
        CC5  
      3.724  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                    Value p-value                Decision
Global Stat        3.3051  0.5081 Assumptions acceptable.
Skewness           0.3324  0.5642 Assumptions acceptable.
Kurtosis           0.0817  0.7750 Assumptions acceptable.
Link Function      1.7360  0.1876 Assumptions acceptable.
Heteroscedasticity 1.1550  0.2825 Assumptions acceptable.

Análisis de varianza

\[Y = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \epsilon_{ijk}\]

\[\hat{Y} = \mu + \tau_{i} + \beta_{j} + \gamma_{k}\]

Para los tratamientos:

\(H_0: \tau_{A} = \tau_{B} = \tau_{C} = \tau_{D} = \tau_{E} = 0\)

\(H_1: \text{En al menos un tratamiento el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un tratamiento.}\)

Para las filas:

\(H_0: \beta_{1} = \beta_{2} = \beta_{3} = \beta_{4} = \beta_{5} = 0\)

\(H_1: \text{En al menos una fila el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_j \neq 0\text{; en al menos una fila.}\)

Para las columnas:

\(H_0: \gamma_{1} = \gamma_{2} = \gamma_{3} = \gamma_{4} = \gamma_{5} = 0\)

\(H_1: \text{En al menos una columna el } \gamma \text{ es diferente a los demás.}\)

\(H_1: \gamma_k \neq 0\text{; en al menos un columna.}\)

anova(modelo.dcl, test = "F")
Analysis of Variance Table

Response: rdt
          Df Sum Sq Mean Sq F value    Pr(>F)    
trt        4 522.30 130.574 27.6685 5.619e-06 ***
F          4  99.20  24.801  5.2553    0.0111 *  
C          4  38.48   9.620  2.0385    0.1527    
Residuals 12  56.63   4.719                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qf(0.95, 4, 12)
[1] 3.259167
shadeDist(qf(0.95, 4, 12), "df",4, 12, lower.tail = F)

Conclusión. A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un tratamiento tiene un efecto sobre el rendimiento estadísticamente diferente del resto de tratamientos.

agricolae::cv.model(modelo.dcl)
[1] 3.984179

Comparaciones de medias

  • A vs B:

\(H_0: \mu_{A} - \mu_{B} = 0\)

\(H_1: \mu_{A} - \mu_{B} \neq 0\)

  • A vs C:

\(H_0: \mu_{A} - \mu_{C} = 0\)

\(H_1: \mu_{A} - \mu_{C} \neq 0\)

  • A vs D:

\(H_0: \mu_{A} - \mu_{D} = 0\)

\(H_1: \mu_{A} - \mu_{D} \neq 0\)

  • A vs E:

\(H_0: \mu_{A} - \mu_{E} = 0\)

\(H_1: \mu_{A} - \mu_{E} \neq 0\)

  • B vs C:

\(H_0: \mu_{B} - \mu_{C} = 0\)

\(H_1: \mu_{B} - \mu_{C} \neq 0\)

  • B vs D:

\(H_0: \mu_{B} - \mu_{D} = 0\)

\(H_1: \mu_{B} - \mu_{D} \neq 0\)

  • B vs E:

\(H_0: \mu_{B} - \mu_{E} = 0\)

\(H_1: \mu_{B} - \mu_{E} \neq 0\)

  • C vs D:

\(H_0: \mu_{C} - \mu_{D} = 0\)

\(H_1: \mu_{C} - \mu_{D} \neq 0\)

  • C vs E:

\(H_0: \mu_{C} - \mu_{E} = 0\)

\(H_1: \mu_{C} - \mu_{E} \neq 0\)

  • D vs E:

\(H_0: \mu_{D} - \mu_{E} = 0\)

\(H_1: \mu_{D} - \mu_{E} \neq 0\)

Prueba de HSD (Honestamente significativa -> Tukey)

agricolae::HSD.test(modelo.dcl, trt = "trt", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

HSD Test for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Alpha: 0.05 ; DF Error: 12 
Critical Value of Studentized Range: 4.50771 

Minimun Significant Difference: 4.379324 

Treatments with the same letter are not significantly different.

     rdt groups
D 59.168      a
E 58.878      a
C 55.728     ab
B 51.718      b
A 47.134      c
agricolae::HSD.test(modelo.dcl, trt = "trt", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

HSD Test for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Alpha: 0.05 ; DF Error: 12 
Critical Value of Studentized Range: 4.50771 

Comparison between treatments means

      difference pvalue signif.        LCL        UCL
A - B     -4.584 0.0388       *  -8.963324 -0.2046756
A - C     -8.594 0.0003     *** -12.973324 -4.2146756
A - D    -12.034 0.0000     *** -16.413324 -7.6546756
A - E    -11.744 0.0000     *** -16.123324 -7.3646756
B - C     -4.010 0.0786       .  -8.389324  0.3693244
B - D     -7.450 0.0012      ** -11.829324 -3.0706756
B - E     -7.160 0.0017      ** -11.539324 -2.7806756
C - D     -3.440 0.1539          -7.819324  0.9393244
C - E     -3.150 0.2126          -7.529324  1.2293244
D - E      0.290 0.9995          -4.089324  4.6693244

Nota: Todos aquellos tratamientos que compartan por lo menos una letra serán estadísticamente similares.

Recomendación: Se debe usar la prueba de Tukey cuando como máximo se tenga 6 tratamientos.

Prueba de Duncan

agricolae::duncan.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = TRUE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

Duncan's new multiple range test
for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Alpha: 0.05 ; DF Error: 12 

Critical Range
       2        3        4        5 
2.993547 3.133384 3.218110 3.274185 

Means with the same letter are not significantly different.

     rdt groups
D 59.168      a
E 58.878      a
C 55.728      b
B 51.718      c
A 47.134      d
agricolae::duncan.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = FALSE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

Duncan's new multiple range test
for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Comparison between treatments means

      difference pvalue signif.        LCL        UCL
A - B     -4.584 0.0059      **  -7.577547 -1.5904531
A - C     -8.594 0.0001     *** -11.727384 -5.4606157
A - D    -12.034 0.0000     *** -15.308185 -8.7598153
A - E    -11.744 0.0000     *** -14.962110 -8.5258900
B - C     -4.010 0.0129       *  -7.003547 -1.0164531
B - D     -7.450 0.0003     *** -10.668110 -4.2318900
B - E     -7.160 0.0003     *** -10.293384 -4.0266157
C - D     -3.440 0.0339       *  -6.573384 -0.3066157
C - E     -3.150 0.0407       *  -6.143547 -0.1564531
D - E      0.290 0.8364          -2.703547  3.2835469

Nota: La prueba de Duncan usa múltiples valores críticos.

Recomendación: No usar la prueba de Duncan cuando se presentan resultados con múltiples tratamientos que obtengan varias significancias.

Prueba de Student - Newman -Keuls

agricolae::SNK.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = TRUE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

Student Newman Keuls Test
for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Alpha: 0.05 ; DF Error: 12 

Critical Range
       2        3        4        5 
2.993547 3.665471 4.079077 4.379324 

Means with the same letter are not significantly different.

     rdt groups
D 59.168      a
E 58.878      a
C 55.728      a
B 51.718      b
A 47.134      c
agricolae::SNK.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = FALSE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

Student Newman Keuls Test
for rdt 

Mean Square Error:  4.719243 

trt,  means

     rdt      std r   Min   Max
A 47.134 2.840129 5 42.26 49.43
B 51.718 5.324328 5 44.41 57.31
C 55.728 2.229747 5 53.72 59.45
D 59.168 2.065798 5 55.87 60.89
E 58.878 1.710181 5 55.87 60.17

Comparison between treatments means

      difference pvalue signif.        LCL        UCL
A - B     -4.584 0.0059      **  -7.577547 -1.5904531
A - C     -8.594 0.0001     *** -12.259471 -4.9285291
A - D    -12.034 0.0000     *** -16.413324 -7.6546756
A - E    -11.744 0.0000     *** -15.823077 -7.6649232
B - C     -4.010 0.0129       *  -7.003547 -1.0164531
B - D     -7.450 0.0008     *** -11.529077 -3.3709232
B - E     -7.160 0.0006     *** -10.825471 -3.4945291
C - D     -3.440 0.0666       .  -7.105471  0.2254709
C - E     -3.150 0.0407       *  -6.143547 -0.1564531
D - E      0.290 0.8364          -2.703547  3.2835469

Prueba de Least significant difference

agricolae::LSD.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = TRUE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

LSD t Test for rdt 

Mean Square Error:  4.719243 

trt,  means and individual ( 95 %) CI

     rdt      std r      LCL      UCL   Min   Max
A 47.134 2.840129 5 45.01724 49.25076 42.26 49.43
B 51.718 5.324328 5 49.60124 53.83476 44.41 57.31
C 55.728 2.229747 5 53.61124 57.84476 53.72 59.45
D 59.168 2.065798 5 57.05124 61.28476 55.87 60.89
E 58.878 1.710181 5 56.76124 60.99476 55.87 60.17

Alpha: 0.05 ; DF Error: 12
Critical Value of t: 2.178813 

least Significant Difference: 2.993547 

Treatments with the same letter are not significantly different.

     rdt groups
D 59.168      a
E 58.878      a
C 55.728      b
B 51.718      c
A 47.134      d
agricolae::LSD.test(modelo.dcl, trt = "trt", 
                       alpha = 0.05, group = FALSE, 
                       main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

LSD t Test for rdt 

Mean Square Error:  4.719243 

trt,  means and individual ( 95 %) CI

     rdt      std r      LCL      UCL   Min   Max
A 47.134 2.840129 5 45.01724 49.25076 42.26 49.43
B 51.718 5.324328 5 49.60124 53.83476 44.41 57.31
C 55.728 2.229747 5 53.61124 57.84476 53.72 59.45
D 59.168 2.065798 5 57.05124 61.28476 55.87 60.89
E 58.878 1.710181 5 56.76124 60.99476 55.87 60.17

Alpha: 0.05 ; DF Error: 12
Critical Value of t: 2.178813 

Comparison between treatments means

      difference pvalue signif.        LCL       UCL
A - B     -4.584 0.0059      **  -7.577547 -1.590453
A - C     -8.594 0.0000     *** -11.587547 -5.600453
A - D    -12.034 0.0000     *** -15.027547 -9.040453
A - E    -11.744 0.0000     *** -14.737547 -8.750453
B - C     -4.010 0.0129       *  -7.003547 -1.016453
B - D     -7.450 0.0002     *** -10.443547 -4.456453
B - E     -7.160 0.0002     *** -10.153547 -4.166453
C - D     -3.440 0.0277       *  -6.433547 -0.446453
C - E     -3.150 0.0407       *  -6.143547 -0.156453
D - E      0.290 0.8364          -2.703547  3.283547

Paquete ExpDes

Diseño completamente al azar

crd.sum <- crd(
  treat = data$trt,
  resp = data$rdt,
  quali = TRUE,
  mcomp = "ccF",
  nl = FALSE,
  hvar = "levene",
  sigT = 0.05,
  sigF = 0.05,
  unfold = NULL
)
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
           DF     SS      MS     Fc      Pr>Fc
Treatament  4 522.30 130.574 13.439 1.7811e-05
Residuals  20 194.32   9.716                  
Total      24 716.61                          
------------------------------------------------------------------------
CV = 5.72 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6338902 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

------------------------------------------------------------------------
Homogeneity of variances test
p-value:  0.1038307 
According to the test of levene at 5% of significance, residuals can be considered homocedastic.
------------------------------------------------------------------------

Test of Calinski & Corsten based on F distribution
------------------------------------------------------------------------
  Groups Treatments  Means
4      a          D 59.168
5      a          E 58.878
3      a          C 55.728
2      b          B 51.718
1      b          A 47.134
------------------------------------------------------------------------

Diseño de bloques completos al azar

dbca.sum <- rbd(
  treat = data$trt,
  resp = data$rdt,
  block = data$`F`,
  quali = TRUE,
  mcomp = "tukey",
  nl = FALSE,
  hvar = "oneillmathews",
  sigT = 0.05,
  sigF = 0.05,
  unfold = NULL
)
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
           DF     SS      MS      Fc     Pr>Fc
Treatament  4 522.30 130.574 21.9656 0.0000025
Block       4  99.20  24.801  4.1721 0.0167511
Residuals  16  95.11   5.944                  
Total      24 716.61                          
------------------------------------------------------------------------
CV = 4.47 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.4781707 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

------------------------------------------------------------------------
Homogeneity of variances test
p-value:  0.02089805 
WARNING: at 5% of significance, residuals can not be considered homocedastic!
------------------------------------------------------------------------

Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    D   59.168 
a    E   58.878 
ab   C   55.728 
 bc      B   51.718 
  c      A   47.134 
------------------------------------------------------------------------

Diseño cuadrado latino

dcl.sum <- latsd(
  treat = data$trt,
  resp = data$rdt,
  row = data$`F`,
  column = data$C,
  quali = TRUE,
  mcomp = "sk",
  sigT = 0.05,
  sigF = 0.05,
  unfold = NULL
)
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
           DF     SS      MS      Fc    Pr>Fc
Treatament  4 522.30 130.574 27.6685 0.000006
Row         4  99.20  24.801  5.2553 0.011101
Column      4  38.48   9.620  2.0385 0.152720
Residuals  12  56.63   4.719                 
Total      24 716.61                         
------------------------------------------------------------------------
CV = 3.98 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6922026 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

Scott-Knott test
------------------------------------------------------------------------
  Groups Treatments  Means
1      a          D 59.168
2      a          E 58.878
3      b          C 55.728
4      c          B 51.718
5      d          A 47.134
------------------------------------------------------------------------

Prueba de Scoott-Knott

Ask <- SK(modelo.dcl,
          which='trt',
          dispersion = 'se',
        sig.level=0.05)
summary(Ask)
Goups of means at sig.level = 0.05 
  Means G1 G2 G3 G4
D 59.17  a         
E 58.88  a         
C 55.73     b      
B 51.72        c   
A 47.13           d

Pruebas no paramétricas

Prueba de Kruskal Wallis

Es aplicable para resultados obtenidos de un DCA.

Requiere del supuesto de que los tratamientos provienen de una misma distribución o tienen distribuciones similares.

data(corn)
str(corn)
'data.frame':   34 obs. of  3 variables:
 $ method     : int  1 1 1 1 1 1 1 1 1 2 ...
 $ observation: int  83 91 94 89 89 96 91 92 90 91 ...
 $ rx         : num  11 23 28.5 17 17 31.5 23 26 19.5 23 ...
psych::histBy(observation~method,
              data=corn)

Se observa que cada método posee una distribución distinta del resto de distribuciones de los métodos.

\(H_0: \text{Las observaciones de los métodos provienen de la misma distribución}\)

\(H_1: \text{Las observaciones de los métodos no provienen de la misma distribución}\)

corn %>% ad.test(observation~method)


 Anderson-Darling k-sample test.

Number of samples:  3
Sample sizes:  34, 34, 34
Number of ties: 64

Mean of  Anderson-Darling  Criterion: 2
Standard deviation of  Anderson-Darling  Criterion: 1.05185

T.AD = ( Anderson-Darling  Criterion - mean)/sigma

Null Hypothesis: All samples come from a common population.

              AD  T.AD  asympt. P-value
version 1: 54.23 49.66        1.150e-30
version 2: 54.60 49.98        5.614e-31

Conclusión: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula. Por lo tanto, las observaciones de los métodos no provienen de la misma distribución.

Recomendación: En base a la evidencia, no se recomienda usar la prueba de Kruskall Wallis.

comparison<-with(corn,
            kruskal(y = observation,
                    trt = method,
                    group=TRUE,
                    main="corn"))
comparison
$statistics
     Chisq Df      p.chisq
  25.62884  3 1.140573e-05

$parameters
            test p.ajusted name.t ntr alpha
  Kruskal-Wallis      none method   4  0.05

$means
  observation     rank      std  r Min Max   Q25  Q50   Q75
1    90.55556 21.83333 3.643869  9  83  96 89.00 91.0 92.00
2    86.40000 15.30000 3.777124 10  81  91 83.25 86.0 89.75
3    95.71429 29.57143 3.638419  7  91 101 93.50 95.0 98.00
4    79.87500  4.81250 1.726888  8  77  82 78.75 80.5 81.00

$comparison
NULL

$groups
  observation groups
3    29.57143      a
1    21.83333      b
2    15.30000      c
4     4.81250      d

attr(,"class")
[1] "group"
comparison<-with(corn,
            kruskal(observation,
                    method,
                    p.adj="bon",
                    group=FALSE, 
                    main="corn"))
comparison
$statistics
     Chisq Df      p.chisq
  25.62884  3 1.140573e-05

$parameters
            test  p.ajusted name.t ntr alpha
  Kruskal-Wallis bonferroni method   4  0.05

$means
  observation     rank      std  r Min Max   Q25  Q50   Q75
1    90.55556 21.83333 3.643869  9  83  96 89.00 91.0 92.00
2    86.40000 15.30000 3.777124 10  81  91 83.25 86.0 89.75
3    95.71429 29.57143 3.638419  7  91 101 93.50 95.0 98.00
4    79.87500  4.81250 1.726888  8  77  82 78.75 80.5 81.00

$comparison
      Difference pvalue Signif.         LCL        UCL
1 - 2   6.533333 0.0426       *   0.1474712 12.9191955
1 - 3  -7.738095 0.0238       * -14.7422174 -0.7339731
1 - 4  17.020833 0.0000     ***  10.2674375 23.7742292
2 - 3 -14.271429 0.0000     *** -21.1206220 -7.4222351
2 - 4  10.487500 0.0006     ***   3.8949224 17.0800776
3 - 4  24.758929 0.0000     ***  17.5658367 31.9520205

$groups
NULL

attr(,"class")
[1] "group"

Prueba de la mediana

out<- with(corn,
           Median.test(observation,
                       method,
                       console=TRUE))

The Median Test for observation ~ method 

Chi Square = 17.54306   DF = 3   P.Value 0.00054637
Median = 89 

  Median  r Min Max   Q25   Q75
1   91.0  9  83  96 89.00 92.00
2   86.0 10  81  91 83.25 89.75
3   95.0  7  91 101 93.50 98.00
4   80.5  8  77  82 78.75 81.00

Post Hoc Analysis

Groups according to probability of treatment differences and alpha level.

Treatments with the same letter are not significantly different.

  observation groups
3        95.0      a
1        91.0      b
2        86.0      b
4        80.5      c
out$statistics
     Chisq Df    p.chisq Median
  17.54306  3 0.00054637     89
out$parameters
    test name.t ntr alpha
  Median method   4  0.05
out$medians
  Median  r Min Max   Q25   Q75
1   91.0  9  83  96 89.00 92.00
2   86.0 10  81  91 83.25 89.75
3   95.0  7  91 101 93.50 98.00
4   80.5  8  77  82 78.75 81.00
out$groups
  observation groups
3        95.0      a
1        91.0      b
2        86.0      b
4        80.5      c
out<- with(corn,
           Median.test(observation,
                       method,
                       console=TRUE,
                       group = F))

The Median Test for observation ~ method 

Chi Square = 17.54306   DF = 3   P.Value 0.00054637
Median = 89 

  Median  r Min Max   Q25   Q75
1   91.0  9  83  96 89.00 92.00
2   86.0 10  81  91 83.25 89.75
3   95.0  7  91 101 93.50 98.00
4   80.5  8  77  82 78.75 81.00

Post Hoc Analysis

        median     chisq pvalue signif.
1 and 2   89.0  2.554444 0.1100        
1 and 3   92.5  6.349206 0.0117       *
1 and 4   83.0 13.432099 0.0002     ***
2 and 3   91.0 13.246753 0.0003     ***
2 and 4   82.5 14.400000 0.0001     ***
3 and 4   82.0 15.000000 0.0001     ***

Prueba de Friedman

data(grass)
attach(grass)
str(grass)
'data.frame':   48 obs. of  3 variables:
 $ judge     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ trt       : Factor w/ 4 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ evaluation: num  4 4 3 3 4 2 1 2 3.5 4 ...
Y <- grass$evaluation
tratamiento <- grass$trt
bloque <- grass$judge

Anp <- friedman(
  judge = bloque,
  trt = tratamiento, 
  evaluation = Y, 
  alpha=0.05, 
  group=TRUE, 
  console=TRUE, 
  main=NULL)

Study: Y ~ bloque + tratamiento 

tratamiento,  Sum of the ranks

      Y  r
t1 38.0 12
t2 23.5 12
t3 24.5 12
t4 34.0 12

Friedman's Test
===============
Adjusted for ties
Critical Value: 8.097345
P.Value Chisq: 0.04404214
F Value: 3.192198
P.Value F: 0.03621547 

Post Hoc Analysis

Alpha: 0.05 ; DF Error: 33
t-Student: 2.034515
LSD: 11.48168 

Treatments with the same letter are not significantly different.

   Sum of ranks groups
t1         38.0      a
t4         34.0     ab
t3         24.5      b
t2         23.5      b